library(tidymodels) #tidymodels 框架
# ✖ purrr::discard() masks scales::discard()
# ✖ dplyr::filter()  masks stats::filter()
# ✖ dplyr::lag()     masks stats::lag()
# ✖ recipes::step()  masks stats::step()
# 执行以上方法需要指定包名运行，如purrr::discard()

library(naniar)   #处理和可视化缺失数据
library(readxl)   #读取EXCEL文件
library(themis)   #处理类别不平衡
library(finetune) #参数调优
library(discrim)  #朴素贝叶斯引擎包
library(themis)   #处理不平衡数据
library(doMC)     #并发处理

# 加载EXCEL数据，数据类型转为numeric
# 读取INR|TG|HDL时存在告警，读取完数据转换应该正常。
# 需要确认为何在EXCEL中的小数无法直接正确预测类型。
mydata = read_excel(paste(getwd(),
                          "data",
                          "aaaalldata.xlsx",
                          sep = "/"),
                    col_types=c("numeric"))
mydata$diagnosis <- as.factor(mydata$diagnosis)

str(mydata) #数据视图，列名|类型|长度|数据
dim(mydata) #数据维度，行数|列数

mydata |> count(diagnosis) # 展示数据比例


# 默认展示10列变量，多余变量通过print(n=...)方式打印
# print(miss_var_summary(mydata),n=42)
miss_var_summary(mydata) #可视化数据缺失比例

# 数据分组，分组方案，7:3，8:2，1:1
set.seed(1)
data_split = initial_split(mydata,
                           prop = 0.7,
                           strata = diagnosis)
train_data = training(data_split)
test_data = testing(data_split)
columns <- c("diagnosis","sex","ALL","N","INR","Cr","bleeding","time","Cervical")
train_data <- train_data[,columns]
test_data <- test_data[,columns]

# 构建预处理方案
base_recp = recipe(diagnosis~.,data = train_data) |>
  step_range(INR,Cr,bleeding,time) |>
  step_smotenc(diagnosis,over_ratio=tune() ,neighbors=tune(),seed=123)


xgboost_spec = boost_tree(tree_depth = tune(),
                          trees = tune(),
                          learn_rate = tune(),
                          mtry=tune(),
                          min_n = tune(),
                          loss_reduction = tune(),
                          sample_size = tune(),
                          stop_iter = tune()) |>
  set_engine("xgboost",counts = FALSE) |>
  set_mode("classification")


base_workflow = workflow() |>
  add_recipe(base_recp) |>
  add_model(xgboost_spec)


base_flow_param = base_workflow |>
  extract_parameter_set_dials() |>
  update(mtry = mtry_prop(c(0.1,1)))

resamples = vfold_cv(data = train_data,
                     v = 5,
                     repeats = 5,
                     strata = diagnosis)

set.seed(123)
registerDoMC(cores = round(parallel::detectCores() * 0.8))
base_tune = tune_bayes(base_workflow,
                       resamples = resamples,
                       initial = 10,
                       iter = 10,
                       param_info = base_flow_param,
                       metrics = metric_set(roc_auc,f_meas,sens,spec),
                       control = control_bayes(verbose = TRUE,
                                               verbose_iter = TRUE,
                                               no_improve = 5,
                                               parallel_over = "everything",
                                               event_level = "second"))

base_tune|>
  show_best(metric = "roc_auc",n=1)
optim_best_auc = base_tune|>
  select_best(metric = "roc_auc")

base_tune|>
  show_best(metric = "f_meas",n=1)
optim_best_f_meas = base_tune|>
  select_best(metric = "f_meas")

base_tune|>
  show_best(metric = "sens",n=1)
optim_best_sens = base_tune|>
  select_best(metric = "sens")

base_tune|>
  show_best(metric = "spec",n=1)
optim_best_spec = base_tune|>
  select_best(metric = "spec")

set.seed(123)
final_auc_workflow  =  finalize_workflow(base_workflow,parameters = optim_best_auc)
best_fit_auc = fit(final_auc_workflow,data=train_data)

set.seed(123)
final_f_meas_workflow  =  finalize_workflow(base_workflow,parameters = optim_best_f_meas)
best_fit_f_meas = fit(final_f_meas_workflow,data=train_data)

set.seed(123)
final_sens_workflow  =  finalize_workflow(base_workflow,parameters = optim_best_sens)
best_fit_sens = fit(final_sens_workflow,data=train_data)

set.seed(123)
final_spec_workflow  =  finalize_workflow(base_workflow,parameters = optim_best_spec)
best_fit_spec = fit(final_spec_workflow,data=train_data)

best_auc_prediction = augment(best_fit_auc,new_data = test_data)
roc_auc(best_auc_prediction,truth = diagnosis,.pred_1,event_level="second")
accuracy(data = best_auc_prediction,truth = diagnosis,.pred_class)
sensitivity(data=best_auc_prediction,truth = diagnosis,.pred_class,event_level = "second")
specificity(data=best_auc_prediction,truth = diagnosis,.pred_class,event_level = "second")
ppv(data=best_auc_prediction,truth = diagnosis,.pred_class,event_level = "second")
conf_mat(data=best_auc_prediction,truth = diagnosis,.pred_class)
f_meas(data=best_auc_prediction,truth = diagnosis,.pred_class,event_level = 'second')


best_f_meas_prediction = augment(best_fit_f_meas,new_data = test_data)
roc_auc(best_f_meas_prediction,truth = diagnosis,.pred_1,event_level="second")
accuracy(data = best_f_meas_prediction,truth = diagnosis,.pred_class)
sensitivity(data=best_f_meas_prediction,truth = diagnosis,.pred_class,event_level = "second")
specificity(data=best_f_meas_prediction,truth = diagnosis,.pred_class,event_level = "second")
ppv(data=best_f_meas_prediction,truth = diagnosis,.pred_class,event_level = "second")
conf_mat(data=best_f_meas_prediction,truth = diagnosis,.pred_class)
f_meas(data=best_f_meas_prediction,truth = diagnosis,.pred_class,event_level = 'second')

best_sens_prediction = augment(best_fit_sens,new_data = test_data)
roc_auc(best_sens_prediction,truth = diagnosis,.pred_1,event_level="second")
accuracy(data = best_sens_prediction,truth = diagnosis,.pred_class)
sensitivity(data=best_sens_prediction,truth = diagnosis,.pred_class,event_level = "second")
specificity(data=best_sens_prediction,truth = diagnosis,.pred_class,event_level = "second")
ppv(data=best_sens_prediction,truth = diagnosis,.pred_class,event_level = "second")
conf_mat(data=best_sens_prediction,truth = diagnosis,.pred_class)
f_meas(data=best_sens_prediction,truth = diagnosis,.pred_class,event_level = 'second')

best_spec_prediction = augment(best_fit_spec,new_data = test_data)
roc_auc(best_spec_prediction,truth = diagnosis,.pred_1,event_level="second")
accuracy(data = best_spec_prediction,truth = diagnosis,.pred_class)
sensitivity(data=best_spec_prediction,truth = diagnosis,.pred_class,event_level = "second")
specificity(data=best_spec_prediction,truth = diagnosis,.pred_class,event_level = "second")
ppv(data=best_spec_prediction,truth = diagnosis,.pred_class,event_level = "second")
conf_mat(data=best_spec_prediction,truth = diagnosis,.pred_class)
f_meas(data=best_spec_prediction,truth = diagnosis,.pred_class,event_level = 'second')



